Dynamical patterns of epidemic outbreaks in complex heterogeneous networks 
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We present a thorough inspection of the dynamical behavior of epidemic phenomena in populations with 
complex and heterogeneous connectivity patterns. We show that the growth of the epidemic prevalence is 
virtually instantaneous in all networks characterized by diverging degree fluctuations, independently of the 
structure of the connectivity correlation functions characterizing the population network. By means of 
analytical and numerical results, we show that the outbreak time evolution follows a precise hierarchical 
dynamics. Once reached the most highly connected hubs, the infection pervades the network in a progressive 
cascade across smaller degree classes. Finally, we show the influence of the initial conditions and the 
relevance of statistical results in single case studies concerning heterogeneous networks. The emerging 
theoretical framework appears of general interest in view of the recently observed abundance of natural 
networks with complex topological features and might provide useful insights for the development of adaptive 
strategies aimed at epidemic containment. 



1. Introduction 

Accurate mathematical models of epidemic spreading 
are the basic conceptual tools in understanding the im- 
pact of diseases and the development of e ffective s trate- 
gies for their control and containment dBajleyl 19751 
Anderson fc Mavl Il992l iDiekmann fc Heesterbeekt hooot 
Dailev fc Ganill200l[L In order to increase the relevance 
and utility of these strategies, it is crucial that models 
contain the key features that appropriately characterize 
the system of interest. The age and social structure of the 
population, the contact network among individuals, and 
the metapopulation characteristics such as the geograph- 
ical patch structure, are all factors that might acquire a 
particular relevance in a reliable epidemic model (for a 
recent review on the subject see ijFerguson et a/.ll2003j) l. 
Among these features, the connectivity pattern of the 
network of contacts among individuals along which the 
disease can be transmitted has been acknowledged since 
a long time as a relevant factor in determinin g the prop- 
erties of epidemic spreading p henom en a (see Yorke^e^^ 
" l978h : iHethcote fc Yorkd lll98l: iMav fc Andersor 
19881) : lAnderson fc Mavl l|l992lk lEubank et all (|2004 
and references therein). In these studies, it was ob- 
served that the heterogeneity of the population network 
in which the disease spreads may have noticeable effects 
in the evolution of the epidemic as well as in the effect of 
immunization strategies. 

The study of the impact of network structure in the 
epidemic modeling has been recently revamped by the 
large number of results on large networked systems point- 



ing out the ubiquitous presence of heterogeneities and 
complex topological features on a wide range of scales 
(Albert & Barabasi, 2002 : iDorogovtsev fc Mendesl 12003: 
IPastor-Satorras fc Vespignanil 12004ft . A striking example 
of this new framework is provided by scale- free networks, 
which are characterized by virtually infinite fluctuations 
in the number of connections k (the degree) that any given 
vertex in the network may have. This feature finds its 
signature in a heavy-tailed degree distribution (defined 
as the probability that any vertex is connected to k other 
vertices) , often approximated by a power-law behavior of 
the form P(k) ~ fc -7 , with 2 < 7 < 3. This implies 
an unexpected statistical abundance of vertices with very 
large degrees; i.e. the s o called "hubs" or "superspread- 
ers" (|Yorke et a l. 1978). In scale-free networks the abun- 
dance of those vertices reaches a level that guarantees 
the proliferation of a large number of infected individu- 
als whatever the rate of infection characterizing the epi- 
demic, eventually leading to the absence of any epidemic 
threshold below which the infection cannot initiate a 
major outbreak llPastor-S atorras k, Vespi gnanil 12001 

a iNewmanl 12 



IMav fc LlovdU200lUMoreno et «dl2002l lNewma.nl l200lh 

This peculiar theoretical scenario, however, turns out 
to be of practical importance since both the sex- 



ual co ntact pattern llColgate et all . Il989t iLilieros et all 
EEl ISchneeherger et all l2004h 



and several techno- 



logical networks llPastor-Satorras fc Vespignanil I200T 



iLlovd fc Mal.l200lUPastor-Satorras fc Vespignani. I20& 
appear to have scale-free features. This implies that 
both sexually transmitted diseases and computer viruses 



1 



may fit in this scenario, raising new questions and 
scrutiny on several epidemic models and strategies 
aimed a t optimizing the deployment of immunization re - 
sources l)Dezso fc BarabasilfeoO^ICohen fc Havlinll2003t) . 

In this paper we provide a thorough presentation 
of results concerning the analysis of the time evolution 
of epidemic outbreaks in complex networks with highly 
heterogeneous connectivity patterns. We consider the 
time behavior of epidemic outbreaks in the general class 
of models without and with internal recovery and find 
that the growth of infected individuals is governed by a 
time scale r proportional to the ratio between the first 
and second moment o f the network's degree distribu- 
tion l|Barthelemv et all Effll . r ~ (k)/(k 2 ). This im- 
plies that the larger the degree fluctuations (governed by 
(fc 2 )), the faster the epidemic propagation will be. In par- 
ticular, a virtually instantaneous rise of the prevalence 
is obtained in scale-free networks where (fc 2 ) — * oo for 
infinite network sizes. This result is shown to be valid 
also in networks with non trivial connectivity correlation 
functions as often encountered in real systems analysis. 
Furthermore, we study the detailed propagation in time 
of the infection through the different degree classes in 
the population. We find a striking hierarchical dynam- 
ics in which the infection propagates via a cascade that 
progresses from higher to lower degree classes. This infec- 
tion hierarchy might be used to develop dynamical ad-hoc 
strategies for network protection. Finally, we study the 
influence of initial conditions on the epidemic develop- 
ment and the relevance of statistical results in the case of 
single case studies. All the analytical discussion is supple- 
mented with careful numerical simulations at the discrete 
individual level. 



2. Basic theory of epidemic dynamics 



2.1. REPRODUCTIVE NUMBER AND TIME EVOLUTION 



where A is the spreading rate, defined as the probability 
rate that a susceptible individual in contact with an in- 
fected individual will contract the disease, and /j, is the 
recovery rate of infected individuals, either to the suscep- 
tible or the recovered states. It is easy to understand that 
any epidemic will spread across a non zero fraction of the 
population only for Rq > 1. In this case the epidemic is 
able to generate a number of infected individuals larger 
than those which are recovered, leading to an increase of 
the infected individuals i(t) at time t following the expo- 
nential form 



i(t) ~ i Q e t/Td . 



(2) 



Here io is the initial density of infected individuals and 
Td is the typical outbrea k time, that in general reads as 
llAnHerson MavLlTflfll 



T^ 1 = h(Rq - 1) 



(3) 



The previous considerations lead to the definition of a cru- 
cial epidemiological concept, namely the epidemic thresh- 
old. Indeed, if the spreading rate is not large enough to 
allow a reproductive number larger than one (A > /i/(fc)), 
the epidemic outbreak will not affect a finite portion of 
the population and dies out in a finite time. In epidemio- 
logical studies for which the two important assumptions of 
"homogeneous mixing" and constant infectiousness (con- 
stant A) are made, the spreading pattern of the epidemics 
is therefore controlled by the generation time scale 1/A 
and i?o an d there are roughly thr ee different stages in an 
epidemics l|Ferguson et all EoOffl. More precisely, when 
infectious individuals are introduced in a network, one 
observes a first noisy phase followed in general by an ex- 
ponential outbreak of the epidemics. Depending on the 
long term behavior of individuals against the disease we 
will observe at large times a different behavior described 
by the specific epidemic model used (see Fig. QJ). 



At first instance epidemiological studies deal with the 
properties of epidemics in the equilibrium or long time 
steady state such as a non-zero prevalence state associ- 
ated to the presence of an endemic phase, the presence or 
absence of a global outbreak, or the non-seasonal cycles 
that are observed in many infections l|Anderson fc Mavl 
1992?) . As well, the dynamical evolution of epidemics out- 
breaks and the effects of the introduction of a seed of in- 
fection into a large population of susceptible individuals 
are of great concern. A basic parameter in epidemiology 
is the basic reproductive number Rq, which counts the 
number of secondary infected cases generated by one pri- 
mary infected individual. Under the assumption of the 
homogeneous mixing of the population if an infected in- 
dividual is in contact with (fc) other individual, the basic 
reproductive number is defined as 

A(fc) 

Ro = — , (i) 

A* 



density of infected 
individuals 




't * s time 

Figure 1 : Typical profile of density of infected versus time on 
a given realization of the network. In the first regime t < ti, 
the outbreak did not start and there are fluctuations. In the 
second regime, ti < t < ti there is an exponential growth 
characterized by Ro. In the final regime (t > t%), the density 
of infected either converges to a constant for the SIS model or 
to zero for the SIR model. 
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The above considerations and parameters are at the 
core of several epidemic models based on the compart- 
mentalization of the population. In other words, each 
individual of the population can only exist in a certain 
number of discrete states such as susceptible, infected or 
permanently recovered. The latter state is equivalent to 
the removal of the individual from the population since 
it is supposed that it cannot get the infection anymore. 
The total population N is assumed to be constant and 
if S(t), I(t) and R(t) are the number of susceptible, in- 
fected and removed individuals at time t, respectively, 
then N = S(t) + I(t) + R(t). The simplest epidemio- 
logical model one can consider is the susceptible-infected- 
susceptible (SIS) model. The SIS model is mainly used as 
a paradigmatic model for the study of infectious diseases 
leading to an endemic state with a stationary and con- 
stant value for the prevalence of infected individuals, i.e. 
the degree to which the infection is widespread in the pop- 
ulation. In the SIS model, individuals can only exist in 
two discrete states, namely, susceptible and infected. The 
disease transmission is described in an effective way. The 
probability that a susceptible vertex acquires the infection 
from any given neighbor in an infinitesimal time interval 
dt is Xdt, where A defines the virus spreading rate. At the 
same time, infected vertices are cured and become again 
susceptible with probability fidt. Individuals thus run 
stochastically through the cycle susceptible — + infected 
— > susceptible, hence the name of the model. The SIS 
model does not take into account the possibility of indi- 
viduals removal due to death or acquired immunization, 
which would lead to t he so-called susceptible - infected- 
remo ved (SIR) model ll Anderson fc Mavt 1 19921: iMurravl 
Il993lh The SIR model, in fact, assumes that infected in- 
dividuals disappear permanently from the network with 
rate \x . In models such as the SIS, the number of in- 
fected individuals increases up to a stationary constant 
value which is non zero if Rq > 1. On the contrary, in 
models such as the SIR, the number of infected individu- 
als tends toward zero since all infected will sooner or later 
become removed from the population. Also in this case, 
however, a finite fraction of the population is affected by 
the epidemic outbreak only if i? > 1 ■ It should be noted 
that it is also possible to induct a steady state in the 
SIR model, by introducing new susceptible individuals at 
a constant rate. This new parameter constitutes a new 
time s cale that gives rise to o scillations in the endemic 
phase l|Mav fc Anderson!, H984j) . 

2.2. COMPLEX HETEROGENEOUS NETWORKS 

The general picture presented in the previous section 
is obtained in the framework of the homogeneous mixing 
hypothesis. This hypothesis assumes that the network of 
contacts among individuals has very small degree fluctu- 
ations. In other words, the degree fc fluctuates very little 



and we can assume k ~ (k), where the brackets (•) denote 
the average over the degree distribution. 

However, networks can be very heterogeneous. So- 
cial heterogeneity and the existence of "super-spreaders" 
have been known for long tim e in the epidemics liter- 
ature l|Hethcote fc Yorkelll984 - The signature of this 
large heterogeneity can be measured in the degree distri- 
bution P(k) which in this case decays very slowly. In- 
deed, in a homogeneou s network such as a random graph 
l)Erdos fe Renvilll959|) . P(k) decays f aster than exponen- 
tially , while for scale-free networks l|Albert fc BarabasiL 
l2002j) it decays as a power law for large k 



P(k) ~ k-' 1 . 



(4) 



Examples of such networ ks relevant to epidemics stud- 
ies in clude the Internet l|Pastor-Satorras fe Vespignani . 
2004 ) , the network o f airline connections ijGuimeraeT^Z , 
200.4 iBarrat et a/1 \2004 or the web of sexual 
contacts llColgate et Qilll989l: iLilieros et all l200ll 
ISchneeberger et all 1200^ . In these networks, the aver- 
age degree (k) is no longer the relevant variable and one 
expects the fluctuations, described by (fc 2 ), to play an 
important role. 

The question of the effect of the heterogene- 
ity on epidemi c behavior has b een a ddressed at 
various levels llHethcote fc YorkeL Il984l) and ana- 
lyzed i n details in the last years f or scale -free net- 
works ijPastor-Satorras fc VespignaniL 12001 tlJ l. These 
studies were concerned with the stationary limit and the 
existence of an endemic phase. A key result is the expres- 
sion of the basic reproductive number which in this case 
takes the form 

(fc 2 ) 



Ro oc 



(fc) 



(5) 



The important fact here is that Rq is proportional to the 
second moment of the degree, which diverges for increas- 
ing network sizes. This has some important epidemiolog- 
ical consequences. Indeed, whatever the spreading rate A 
the basic reproductive rate is always larger than one, thus 
leading to the lack on any epidemic threshold. In other 
words, in heterogeneous networks, whatever the infection 
rate, the epidemics has a finite probability to generate a 
major outbreak. 

These results, while extremely important, consider 
only one face of networks complexity. In general, how- 
ever, real networks have other complex features such as 
strong degree correlation among connected vertices. This 
feature is mathematically characterized through the con- 
ditional probability P(fc'|fc) that a vertex of degree fc is 
connected to a vertex of degree fc'. Networks which are 
completely defined by the degree distribution P(k) and 
the condi tional probability -P(fc'lfc) are call ed Markovian 
networks l|Boguna fc Pastor-Sato7rasll2002j) and must ful- 
fill the following degree detailed balance condition 



kP{k'\k)P{k) = k'P(k\k')P(k'). 



(0) 
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This expression is a mathematical statement of the obvi- 
ous observation that, in any real network, all edges must 
point from one vertex to another. 

The full knowledge of the function P(k'\k), which 
measures the correlations in the network, would of- 
ten be difficult to interpret. Therefore, the aver- 
age nearest neighbor degree k nn , and the behavior of 
this quantity as a function of the degree, k nn (k) — 
J2 k , k'P( k'\k), have been proposed to measure these cor- 
relations l|Pastor-Satorras fc Vespjgnanil . 12004b In the 
absence of degree correlations, P(k'\k) does not depend 
on k and neither does the average nearest neighbors' de- 
gree; i.e. k nn (k) = const. llPastor-Satorras k. VespignaniL 
2004J). In the presence of correlations, the behavior of 
k n n{k) identifies two general classes of networks. If k nn {k) 
is an increasing function of fc, vertices with high degree 
have a larger probability to be connected with large degree 
vertices. This property is refe rred in physics a nd social 
sciences as assortative mixing llNewmanl Enni. On the 
contrary, a decreasing behavior of k nn {k) defines disas- 
sortative mixing, in the sense that high degree vertices 
have a majority of neighbors with low degree, while the 
opposite holds for low degree vertices. These possibilities 
are summarized in Figure [3 Among real networks, many 
social networks display assortative mixing, while techno- 
logical networks show typically disassortative properties. 
In summary, the average nearest neighbor degree carry an 
information about two-point correlations which is easy to 
interpret and avoids the fine details of the full distribution 
P{k'\k). 



k J k > 




sassortative 



Figure 2: Typical behavior of the average nearest neighbor 
degree k nn (k) for assortative mixing and disassortative mix- 
ing. 



The propagation of epidemics on a network may a 
priori be modified by the presence of correlations. The 
quantity k nn seems in this context very relevant, since it 
measures the number of individuals that may be reached 
by the infection in few steps. In gen eral, however, it has 
been shown in l|Boguna et al 1200.^ that, in any scale- 
free network with diverging second moment (fc 2 ), the epi- 
demic threshold of the model vanishes, regardless of the 
presence of correlations. In finite networks, the role of 



correlations appears through the fact that the epidemic 
threshold is b ounded from above b y a quantity which de- 
pends on k nn i Boguna et al 12008^ . 

While these results contribute to a general under- 
standing of epidemic spreading in complex networks, they 
do not provide a large amount of information about the 
time patterns of the dynamical process and how effec- 
tively the spreading occurs. Uncovering these patterns 
and their relation to the heterogeneity of the networks is 
crucial for assessing control strategies. In the following 
sections we will approach this problem within a general 
formalism, valid for most kind of epidemic processes in 
heterogeneous networks. 

3. The early stages of the epidemic outbreak 

In order to understand how the topology of the net- 
work of contacts affects the dynamical process of epi- 
demics spreading, we will perform a study on different 
kinds of networks. Most complex networks can be clas- 
sified according to the d ecay of their degree distribu- 
tion ijAmaral et all l2000|) . We will first consider the 
case where the degree distribution decays faster than any 
power-law, typically as an expone ntial (or faster), such a s 
for the random (Poisson) graph (|Erdos fc RenviL Il959l) . 
For this sort of networks, the degree fluctuations are very 
small and we can approximate the degree of every ver- 
tex as a constant, k ss (k). This approximation used in 
most previous studies l|Anderson fc Mavl Il992l : iMurravi 
corresponds to the homogeneous mixing case in 
which all individuals have the same environment, and the 
same number of acquaintances that can be infected. We 
will then consider the case of heterogeneous networks, ex- 
emplified in the case of the scale-free networks, and show 
how the large heterogeneity induces changes in the epi- 
demic process. 

Both the SIS and SIR models introduce a time scale 
r r = I//i governing the self-recovery of individuals. Since 
we are especially interested in the onset dynamics of out- 
breaks two scenarios are in order. If r r is smaller than 
the spreading time scale, then the process is dominated 
by the natural recovery of infected to susceptible or re- 
moved individuals. This situation is however less inter- 
esting since it corresponds to a dynamic process gov- 
erned by the decay into a healthy state and the inter- 
action with neighbors plays a minor role. Epidemiolog- 
ical concern is therefore in the regime r r 3> 1/A; i.e. a 
spreading time scale much smaller than the recovery time 
scale. In this case, as a first approximation, we can ne- 
glect the individual recovery that will occur at a much 
later stage and focus on the early dynamics of the epi- 
demic outbreak. This case corresponds to the simpli- 
fied susceptible-infected (SI) model, for which infected 
nodes remain always infective and spread the infection 
to susceptible neighbors with rate A; the seed of infec- 
tives placed at time t = will thus infect the rest of 
the network, the dynamical process being controlled by 
the topology of the network. The SI model thus appears 
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as the simplest framework for assessing the effect of net- 
work's topology on the spreading dynamics. This argu- 
ment is the reason why we will present our calculations 
for the simple SI model. For the sake of completeness, 
however, we will readily show in the next sections that 
the results obtained for the SI model can be considered 
as fairly general and can be easily extended to SIS and 
SIR models for heterogeneous networks. 

3.1. HOMOGENEOUS NETWORKS 

A first analytical description of the SI model can 
be u ndertaken within t he homogeneous mix ing hypoth- 
esis l|Anderson fc Mavi Il992t iMurravl Il993h . consisting 
in a mean-field description of the system in which all ver- 
tices are considered as being equivalent. In this case the 
system is completely defined by the number of infected 
individual I(t), and the reaction rate equation for the 
density of infected individuals i(t) = I(t)/N (where N is 
the total size of the population) reads as 



di(t) 
dt 



= X(k)i(t) [1 - i(t)} 



(7) 



The above equation states that the growth rate of infected 
individuals is proportional to the spreading rate A, the 
density of susceptible vertices that may become infected, 
s(t) = 1 — i(t), and the number of infected individuals in 
contact with any susceptible vertex. The homogeneous 
mixing hypothesis considers that this last term is sim- 
ply the product of the number of neighbors (k) and the 
average density i(t). Obviously, this approximation ne- 
glects correlations among individuals and considers that 
all vertices have the same number of neighbors (fc), i.e. it 
assumes a perfectly homogeneous network. The solution 
of equation J2J reads 



i(t) 



i exp(t/r H ) 



1 + i [exp(t/T H ) - 1] 



(8) 



where io is the initial density of infected individuals and 
th = (A(fc)) -1 is the time scale of the infection growth. 
At small times, when the density of infected vertices is 
very small, the leading behavior is given by i(t) ~ ioe*/ Tff , 
i.e. by an exponential growth. The value of th corre- 
sponds to the intuitive fact that small spreading rates 
decrease the growth velocity, while a larger number of 
neighbors ensures a faster invasion of the network. 

3.2. HETEROGENEOUS NETWORKS 

The above calculations are valid for networks in which 
the degree fluctuations are very small, i.e. k « (k) for 
all vertices in the network. We face however a different 
situation in networks with a heterogeneous connectivity 
pattern. In this case the degree k of vertices is highly fluc- 
tuating and the average degree is not anymore a mean- 
ingful characterization of the network properties. In or- 
der to take this into account, it is possible to write down 



the reaction rate equations for the densities of infected 
vertices of degree k, i k (t) = I k (t)/N k , where N k (resp. 
Ife(t)) is the number of v ertices (resp. infected vertices) 
within each degree class k l|Pastor-Satorras & VesnignaniL 
12001^ . In the case of the SI m odel the evolution equa- 
tions read l)Boguna et all l2003|) 



dt 



= X[l-i k (t)} kQ k (t), 



(9) 



where the creation term is proportional to the spreading 
rate A, the degree k, the probability 1 — i k that a vertex 
with degree k is not infected, and the density of in- 
fected neighbours of vertices of degree k. The latter term 
is thus the average probability that any given neighbor of 
a vertex of degree k is infected. 

The simplest situation one can encounter corresponds 
to a complete lack of correlations (We will address the 
case of correlated networks in section 3.4). A network is 
said to have no degree correlations when the probability 
that an edge departing from a vertex of degree k arrives 
at a vertex of degree k' is independent of the degree of 
the initial vertex k. In this case, the probability that each 
edge of a susceptible is pointing to an infected vertex of 
degree k' is proportional to the fraction of edges emanated 
from these vertices. By considering that at least one of 
the edges of each infected vertex is pointing to another 
infected vertex, from which the infection has been trans- 
mitted, one obtains 



Q k (t) = Q(t) 



EAk'-l)P(k')i k ,(t) 
(k) 



(10) 



where (k) = ^2 k , k! P(k') is the proper normalization fac- 
tor dictated by the total number of edges. At t = 0, if 
the initial density is not too large, all the neighbors of the 
initially infected vertices are susceptible and we have 



e(t = o) = 



Eg k'P(k<)i k ,(0) 

(k) 



(ii) 



A reaction rate equation for O(i) can be obtained from 
Eqs. © and (jTU|) . In the initial epidemic stages, neglect- 
ing terms of order 0(i 2 ), the following set of equations is 
obtained 



di k (t) 

dt 
dG(t) 

dt 



= XkQ(t), 

(k) 



A 



i e(t). 



(12) 
(13) 



These equations can be solved, yielding for the prevalence 
of nodes of degree k, in the case of an uniform initial con- 
dition i k (t = 0) = io, 



i k {t) = io 



1 



(k 2 ) ~ (k) 



(14) 



and for the total average prevalence i(t) = J2 k P{k)i k (t) 

\2 



i(t) = i 



1 



(k 2 )-(k) [G ij 



(15) 
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where 



(k) 



' \{(k*)-{k))- (16) 

The result Eq. (JTSJ) for uncorrelated networks readily 
implies that the growth time scale of an epidemic out- 
break is related to the graph heterogeneity. Indeed, the 
ratio 



(k) 



(17) 



is the parameter defining the level of heterogeneity of the 
network, since the normalized degree variance can be ex- 
pressed as n/{k) — 1 and therefore high levels of fluctu- 
ations correspond to k > 1. In homogeneous networks 
with a Poisson degree distribution, in which k = (k) + 1 
we recover the result r = (A(fc)) -1 , corresponding to the 
homogeneous mixing hypothesis. On the other hand, in 
networks with very heterogeneous connectivity patterns, 
k is very large and the outbreak time-scale r is very small, 
signaling a very fast diffusion of the infection. In partic- 
ular, in scale-free networks characterized by a degree ex- 
ponent 2 < 7 < 3 we have that k ~ (fc 2 ) — ► oo with the 
network size N — > oo. Therefore in uncorrelated scale- 
free networks we face a virtually instantaneous rise of the 
epidemic incidence. 

3.3. EXTENSION TO THE SIS AND SIR MODELS 

It is worth stressing that the above results can 
be easily extended to the SIS and the SIR models 
l|Anderson fc Mavl [i992h In the case of uncorrelated net- 
works, Eq. © contains, for both the SIS and the SIR 
models, an extra term —/j,i k (t) defining the rate at which 
infected individuals of degree k recover and become again 
susceptible or permanently immune and thus removed 
from the population, respectively: 



dik(t) 
dt 



= Xks k (t)Q k (t) - frik(t) 



(18) 



In the SIS model we have, as usual, s k (t) = 1 — i k (t). 
In the SIR model, on the other hand, the normalization 
imposes that s k (t) = 1 — i k (t) — r k (t), where rk(t) is the 
density of removed individuals of degree k. The inclusion 
of the decaying term —/j/i k , does not change the picture 
obtained in the SI model. By using the same approxima- 
tions, the time scale is found to behave as 



(k) 



A(fc 2 ) - {n + X){k) 



(19) 



In the case of diverging fluctuations the time-scale behav- 
ior is therefore still dominated by (fc 2 ) and r is always pos- 
itive and going to as (fc 2 ) — > oo whatever the spreading 
rate A. This allows to recover the absence of an epidemic 
threshold, i.e. the lack of a decreasing prevalence region 
in the parameter space. It must be noted, however, that 
if fluctuations are not diverging, the outbreak time scale 
is slightly different for the SIS and SIR models (equation 
mi} ) than for the SI model (equation lflE)l ). 



3.4. HETEROGENEOUS CORRELATED NETWORKS 

While we have so far restricted our study 
to the case of uncorrelated networks, it is 
worth noting that many real networks d o not 
fill this assumption foorogovtsev fc Mendesl |2003; 
IPastor-Satorras fc VespignaniL 12004^ In order to con- 
sider the presence of non trivial correlations we have to 
fully take into account the structure of the conditional 
correlation function P(fc'|fc). The equations we have writ- 
ten for the e volution oi i k in the SI model can therefore 
be stated as i|Boguna et all l2003|) 



dt 



x[i-i k (t)]ke k (t) 
fc' - i 



k' 



-P{k'\k). 



(20) 



Here the Qk function takes into account explicitly the 
structure of the conditional probability that an infected 
vertex with degree fc' points to a vertex of degree fc, with 
any of the fc' — 1 free edges it has (not pointing to the orig- 
inal source of its infection). In the absence of correlations, 
it is possible to see that P(k'\k) — fc'P(fc')/(fc), recovering 
the results of section 3.2. If the network presents correla- 
tions, measured by P(k'\k), the situation is slightly more 
complex. The rate equation for i k (t) can be written in 
this case, neglecting terms of order 0(i 2 ), as 

dik[t) = j:xk^p(k'\k)i k ,(t) 



dt 



Ck.k'ik' (t), 



(21) 



which is a linear systems of differential equations given 
by the matrix C = {C k ,k'} of elements 



C k ,k> =Afc^— -P{k'\k) 



(22) 



Elementary considerations from mathematical analysis 
tell us that the behavior of i k (t) will be given by a 
linear combination of exponential functions of the form 
exp(Ajt), where A< are the eigenvalues of the matrix C. 
Therefore, the dominant behavior of the averaged preva- 
lence will be 



iit) 



„A m t 



(23) 



where A m is the largest eigenvalue of the matrix C. 
In the case of an uncorrelated network, we have that 
C k lc k , — Afc(fc' — l)P(fc')/(fc), which has a unique eigen- 
value satisfying 



^C fc , fe ^ = A^* 



(24) 



where = A((fc 2 )/(fc) — 1), and where the correspond- 
ing eigenvector is ^k — k, thus recovering the previous 
result Eq. I|ltj[l for this kind of networks. 
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In the case of correlated networks, it has be en shown 
using the Frobenius theorem l|Gantmacheiill97^ that the 
largest eigenvalue is bounded from below l|Boguha et al 
2003) 

A^>min^^(fc'-l)G-l)P(/|fc)P(fc'|0 . (25) 

k' I 

This equation is very interesting since it can be rewritten 
as 

A 2 m > min]T(Z - l)P(l\k)(k nn (l) ~ 1) ■ (26) 
i 

It has been shown ijBoeuna et al l2003j) that, for scale- 
free network with 2 < 7 < 3, k nn (l) diverges for infinite 
size systems (N —> 00). This ensures that also A TO di- 
verges. Two particular cases have however to be treated 
separately: it may happen that, for some fc , P(l\ko) = 0; 
then the previous limit for gives no information but it 
is then possible to show with slightly more involved calcu- 
lation s that A m still diverges for N — > oo (|Boguha et oil 
120031) . Another problem arises if k nn {l) diverges only for 
1 = 1; this happens however only in particular networks 
where the singularity is accumulated in a pathological way 
onto vertices with a single edge. Explicit examples o f this 
situation are provided in l|Vazauez fc Morenoll2003|) . 

The previous result has the important consequence 
that, even in the presence of correlations, the time scale 
t ~ 1/Am tends to zero in the thermodynamic limit for 
any scale- free network with 2 < 7 < 3. It also underlines 
the relevance of the quantity k nn , which gives a lower- 
bound for A m in finite networks. In the next sections we 
will analyze numerically both the SI and SIS models in 
order to provide a full account of the dynamical proper- 
ties that takes into account the network's complexity as 
well as finite size effects in the population. 

4. Numerical simulations 

In order to test the analytical predictions made in the 
previous sections, we have performed extensive numerical 
simulations of the SI and the SIS model in two different 
paradigmatic examples of complex network models with 
homogeneous and heterogeneous properties. The choice of 
the network models was dictated by the request of gener- 
ating an uncorrelated network, to fully exploit the analyt- 
ical predictions. Simulations use an agent-based modeling 
strategy in which at each time step the SI dynamics is ap- 
plied to each vertex by considering the actual state of the 
vertex and its neighbors. It is then possible to measure 
the evolution of the number of infected individuals and 
other quantities. In addition, given the stochastic nature 
of the model, different initial conditions and networks re- 
alizations can be used to obtain averaged quantities. In 
our simulations, we use N = 10 4 and (k) ranging from 4 
to 20. We typically average our results for a few hundred 
networks while for each network, we average over a few 
hundreds different initial conditions. 



4.1. HOMOGENEOUS NETWORKS 

The example of homogeneous complex network we 
have chosen is t he random graph model proposed b y 
Erdos and Renyi l)Erdos k Renvii Il959t iBollobast Il985|) . 
The network is constructed from a set of N different ver- 
tices, in which each one of the N(N — l)/2 possible edges 
is present with probability p (the connection probability) , 
and absent with probability 1 — p. This procedure results 
in a random network with average degree (k) = pN and 
a Poisson degree distribution in the limit of large N and 
constant (k), 

P(k) = e-^^. (27) 

For (k) > 1, it can be proved that the network exhibits a 
giant component, that is, a set of connected vertices whose 
size is proportional to N. In an actual realization, how- 
ever, it is possible to generate networks in which a fraction 
of the vertices belong to disconnected clusters. Therefore, 
we make the computer simulations of the spreading mod- 
els only on the giant component. 

In the case of the random graph, we checked the va- 
lidity of equation (JHJ , and that the time scale is given by 
1/X{k), as can be seen in Fig. ©. 




Figure 3: Main frame: the symbols correspond to simula- 
tions of the SI model with A = 10 -4 on ER networks with 
N = 10 4 , (k) = 20, 40; the lines are fits of the form of Eq. ©. 
Inset: measured time scale t, as obtained from fitting, versus 
the theoretical prediction for different values of (k) and A. 



Needless to say, in the case of a homogeneous network, 
the hypothesis k ~ (k) captures the correct dynamical be- 
havior of the spreading. This is a standard result and we 
report the numerical simulations just as a reference for 
comparison with the following numerical experiments on 
heterogeneous networks. 

4.2. HETEROGENEOUS NETWORKS 

As a typical example of heterogeneous network, we 
have chosen the network s generated with the Bar abasi- 
Albert (BA) algorithm (|Barabasi fc Albert! Il999lh In 
this algorithm, the network starts from a small core of 
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mo connected vertices. At each time step a new ver- 
tex is added, with m edges (m < mo) connected to the 
old vertices in the network. The new edges are con- 
nected to an old vertex i with a probability proportional 
to its degree fcj (preferential attachment). The networks 
generated by this algorithm have a minimum degree m, 
an average degree (k) = 2m, a scale-free degree distri- 
bution P(k) ~ fc~ 3 and vanishingly small correlations 
l|Pastor-Satorras fc VesDignanil Eo04'l. In our simulations 
we use different network sizes N and minimum degree val- 
ues m in order to change the level of heterogeneity, that 
in this case is given by k ~ mlniV. 

In the simulation we keep track of the average density 
of infected individuals versus time for a network generated 
with the preferential attachment rule. In Fig. 0] we show 
a typical result and the fitting procedure which allows us 
to measure the value of r. 

lr 




Figure 5: Measured time scale r in BA networks as obtained 
from exponential fitting versus the theoretical prediction for 
different values of m and N corresponding to different levels 
of heterogeneity. 

4.3. THE SIS DYNAMICS 

As claimed in previous sections, we have tested that at 
short times and for large heterogeneous networks (k» 1), 
the SI model is a good approximation to the more general 
SIS model described by Eq. fl%|l . In Fig. we plot i(t) 
versus t for a simulation of the SI and the SIS model with 
A = 10~ 3 , on a BA network with m = 8 and N = 10 4 , 
corresponding to (k) = 16 and (k 2 ) sw 632. We took 
different values for the parameter \x: fi — (which corre- 
sponds to the SI model), fi = 1(T 4 , 1(T 3 , 1CT 2 ; Eq. HU 
yields the following values for r: r(/x = 0) w 26, t(/i = 
0.0001) « 26.1, t(jj, = 0.001) w 26.7, r(fi = 0.01) w 35. 
Figure clearly shows that indeed for fi(k) <C A(fc 2 ) and 
short times the curves are indistinguishable so that the 
SIS model is well described by the SI approximation in 
this regime. As fi gets larger the behavior changes since 
the time scale of the recovery is no longer much larger 
than the other time scales. It is possible to show that 
that the same results apply also in the case of the SIR 
model. 



1000 2000 3000 4000 5000 6000 

t 



Figure 4: Average density of infected individuals versus time 
in a BA network of N=10 4 with m = 2. The inset shows the 
exponential fit obtained in the early times (lines) and the nu- 
merical curves i(t) for networks with m = 4, 8, 12, 20 (from 
bottom to top). 

In Fig. we report the early time behavior of out- 
breaks in networks with different heterogeneity levels 
and the behavior of the measured r with respect to 
(k)/X((k 2 ) — (k)). The numerical results recover the an- 
alytical prediction with great accuracy. Indeed, the BA 
network is a good example of uncorrelated heterogeneous 
network in which the approximations used in the calcu- 
lations are satisfied. In networks with correlations we 
expect to find different quantitative results but a qualita- 
tively similar framework as it happens in the case of the 
epidemic threshold evaluation l|Boguha et adl2003Th 





2000 



Figure 6: SIS model on a BA network with N = 10 4 , m = 8; 
A = 1CT 3 , and u = 0, 1CT 4 , 1CT 3 , 1(T 2 from top to bottom. 
The inset focuses on the short time dynamics: the curve 
is superimposed on data for the SI model up to more than 
10r(/x = 0) for fi — 10 -4 and up to a few times r(fi = 0) for 
u = 1(T 3 . 



<k>/X(<k > - <k>) 



5. The infection time pattern: the cascade effect 

The previous results show that the heterogeneity of 
scale-free connectivity patterns favors epidemic spread- 
ing not only by suppressing the epidemic threshold, but 
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also by accelerating the epidemic propagation in the pop- 
ulation. The velocity of the spreading leaves us with 
very short response times in the deployment of control 
measures and a detailed knowledge of the way epidemics 
propagate through the network could be very valuable in 
the definition of adaptive strategics. Indeed, the epidemic 
diffusion is far from homogeneous. The simple formal in- 
tegration of Eq. © written for Sk — 1 — %k yields 



Sk(t) 



= s°e" Afe<I>(t) 



(28) 



where = J Q dt'Q(t'). This result is valid for any 

value of the degree k and the function $(£) is positive 
and monotonously increasing. This last fact implies that 
Sk is decreasing monotonously towards zero when time 
grows. Thus, if one has two values k > k 1 and whatever 
the initial conditions s° and s^,, there is a time t x after 
which Sk{t) < Sk'(t). This 'crossing' time is given by 



0(*x) - 



1 



A(fc - k') 



log 



(29) 



which means that the vertices with the lowest degree are 
typically the last to be infected. 




Figure 7: Time behavior of the average degree of the newly 
infected nodes for SI outbreaks in BA networks (here of size 
N — 10 4 ). Time is rescaled by r. Reference lines are drawn 
at the asymptotic values {k 2 } / {k) for t <^ r and m for Or. 
The two curves are for m = 4 (bottom) and m = 14 (top). 



(If Sj! < SjL, there is no crossing and for all times 
Sfc(t) < Sk'(t)). This result indicates that after an ini- 
tial regime which depends on the initial conditions, the 
disease spreads always from large connectivity to smaller 
connectivities. 

We can go beyond the analytical result that is neces- 
sarily grounded on average quantities by relying on nu- 
merical simulations. Indeed, a more precise characteri- 
zation of the epidemic diffusion through the network can 
be achieved by studying some convenient quantities that 
highlight the invasion pattern of the infection in numerical 
spreading experiments in BA networks. First, we measure 
the average degree of the newly infected nodes at time t, 
defined as 



fcinf(*) 



I(t)-I(t-1) 



(30) 



In Fig. [7] we plot this quantity for BA networks as a func- 
tion of the rescaled time t/r. The curves show an initial 
plateau that can be easily understood by considering that 
at very low density of infected individuals, each vertex 
will infect a fraction of its neighbors without correlations 
with the spreading from other vertices. In this case each 
edge points to a vertex with degree k with probability 
kP{k)/(k) and the average degree of newly infected ver- 
tices is given by 



AtoW = (k 2 )/(k) 



(31) 



After this initial regime, ki n { (t) decreases smoothly when 
time increases. The dynamical spreading process is there- 
fore clear: after the hubs are very quickly infected the 
spread is going always towards smaller values of k. This 
is confirmed by the large time regime that settles in a 
plateau 

feinf (t) = m , (32) 



Further information on the infection propagation 
is provided by the i nverse participation rat io Y?(t ) 



l)Perrida fc Flvvbieret Il987t iBarthelemv et all EoO 



We first define the weight of infected individuals in each 
degree class k by Wk(t) = Ik(t)/I(t). The quantity Y% is 
then defined as 



k 



w 2 k (t). 



(33) 



If Y% ~ ^/k ma x (kmax is the maximal connectivity), in- 
fected vertices are homogeneously distributed among all 
degree classes. In contrast, if Y% is not small (of order 
l/n with n of order unity) then the infection is localized 
on some specific degree classes that dominate the sum 
of Eq. (|33|l . In Fig. [S] we report the behavior of Y 2 ver- 
sus time for BA networks with different minimum degree. 
The function Y% has a maximum at the early time stage, 
indicating that the infection is localized on the large k 
classes, as we infer from the plot of fci n f(i), see Fig. [7| 
First Y 2 decreases, with the infection progressively invad- 
ing lower degree classes, and providing a more homoge- 
neous diffusion of infected vertices in the various k classes. 
Finally, the last stage of the process corresponds to the 
capillary invasion of the lowest degree classes which have a 
larger number of vertices and thus provide a larger weight. 
In the very large time limit, when the whole network is 
infected, Y 2 {t — 00) = J2k^(k) 2 - Noticeably, curves for 
different levels of heterogeneity have the same time pro- 
file in the rescaled variable t/r. This implies that, despite 
the various approximations used in the calculations, the 
whole spreading process is dominated by the time-scale 
defined in the early exponential regime of the outbreak. 

The presented results provide a clear picture of the 
infection propagation in heterogeneous networks. First 
the infection takes control of the large degree vertices in 
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the network. Then it rapidly invades the network via a 
cascade through progressively smaller degree classes. The 
dynamical structure of the spreading is therefore charac- 
terized by a hierarchical cascade from hubs to intermedi- 
ate fc and finally to small fc classes. 




Figure 8: Inverse participation ratio Yi versus time for 
BA network of size N = 10 4 with minimum degree m = 
4, 6, 8, 10, 12, 14 and 20, from top to bottom. Time is rescaled 
with r . The reference line indicates the minimum of Yi around 
t/r ~ 6.5. 



6. The effect of the initial seed of the infection 

In very heterogeneous networks it is reasonable to 
forecast that epidemic outbreaks starting from individ- 
uals possessing very different connectivity properties may 
undergo a rather different time evolution. In order to 
investigate in more details the effect of initial condi- 
tions, one can write the general solution of Eqs. (|13|) with 
ik{t = 0) = /(fc) where / is some given function of the 
degree 



**(*) = f(k) + kG 



(k) 



(fc 2 > - <fc> 



1 



(34) 



where Oo = 6(£ = 0) = (kf(k))/(k). Comparing this last 
expression with Eq. 1)14(1. we conclude that, apart from 
the effect of the initial density of infected individuals, the 
only effect in the dynamics is the presence of the prefactor 
0o, which implies that the spreading is faster for larger 
6o- In the particular case in which the initial infection 
is located on No vertices of given degree fco, we obtain 
i ko (t = 0) = N /(NP(k )) and thus 6 - k N Q /(N(k)). 
Analogously it is easy to see that, for A*o initially contami- 
nated sites, all with connectivity larger than a certain fcjw, 
0o is proportional both to Nq/N and to fcjvf. 

Eq. (|34|) clearly highlights the role of the large con- 
nectivity sites: on the one hand, for a given initial condi- 
tion, the infection of classes with larger k leads to faster 
rise of the epidemic; on the other hand, the larger the 
degree of the initially infected sites, the faster the propa- 
gation of the epidemics. Note however that the time-scale 



r itself is not affected by the initial conditions, which ap- 
pear through the prefactors only. Using Eq. (|34|) we can 
highlights these features since the right hand side of the 
expression 



ik(t) ~ f(k) 
kS 



(k) 



(k 2 ) - (k) 



(35) 



has to be independent on both the initial conditions 
(which enter only through f(k) and Oq) and the degree 
k. We have verified this in Fig. by means of numerical 
simulations of the SI model on BA networks with specific 
initial conditions: if the spreading starts from Nq infected 
sites of connectivity fco, Qo is proportional to Nq and fco, 
so that we show that the quantity 



ik(t)-f(k) 
kN k 

is independent of fc and fco, and is also equal to 

i(t) - io 
(k)N k ' 



(36) 



(37) 



At large times the curves separate since the expression 
<|14|) is no longer valid. The data collapse of Fig. [5] once 
again shows the validity of the analytical approach devel- 
oped in section 3. The time scale r for the exponentially 
fast outbreak of the epidemics is relevant for any kind of 
initial conditions and for all nodes, independently from 
their connectivity. The prefactors however show that the 
spreading is faster if the infected seeds have a large con- 
nectivity and, among the connectivity classes, hubs are 
typically infected faster. 




Figure 9: Collapse plot of quantities l|3lijl and IpTf)) vs. 
time for the SI model on a BA network, with initial con- 
ditions given by No infected sites having a given connec- 
tivity fc - Here N Q = 10, N = 10 4 . The data are av- 
eraged over 200 realizations of the network and 50 runs 
on each network sample. Lines correspond to 1)36(1 for 
fco = 10, 20, 30, 50; symbols correspond to lj3*7|) for the 
same initial conditions and various values of fc between 10 
and 90 (The maximal connectivity is of order y/~N = 100). 

7. Fluctuations: the relevance of t for single case 
studies 
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The previous results are valid for the average density 
of infected nodes and it is legitimate to ask about the rel- 
evance of these results to the case of a single network sam- 
ple. Indeed, in the real-world there is not such a thing as 
averages over different network realizations and we have 
to check the robustness of our results on a single network 
sample. It is then natural to wonder if large statistical 
fluctuations in the network structure may reverberate on 
the difference between single case studies and asymptotic 
average results. In particular, we can compare the pre- 
viously defined time scale r (Eq. which corresponds 
to the exponential growth of the prevalence averaged over 
the networks statistical properties and initial conditions, 
with an operative measure of time-scale obtained in a sin- 
gle network realization case study. For a given network, 
and starting from a given initial density of infected ver- 
tices «o, we define for each run of the spreading process 
the typical realization time-scale as the time r' at which 
the density of infected nodes i run (t) in the specific run is 
equal to the average quantity i(t) at time r provided by 
Eq. I|15|) . This defines r' as the effective time at which 
i r un{t = t') = i{t = t) (see Figure ITU)) . The deviation 
from the average behavior in each run can be readily ac- 
counted by the variable u — t' /t which is a simple mea- 
sure of the relevance of r for each run. In Fig. Illf a - ). 
we plot the probability distribution of u computed for a 
given realization of the network. This figure shows that 
P(u) is not broadly distributed. In other words it is very 
unlikely to find large deviations of u from the unity, at 
which the distribution is strongly peaked. This evidence 
indicates that even if fluctuations are possible, the time- 
scale defined by r is a meaningful characterization of the 
spreading process. 



density of infected 
individuals 




time 



Figure 10: Schematic definition of r' for a single run. Given 
the average density i(t), r' is defined as the time in which the 
density of a given run reaches i(r). 



A more detailed analysis highlights once more the key 
role of initial conditions in the epidemics spread. In 
Fig- HTT b). we used the same network as in Fig. ITlT a^. 
but we made the following distinction according to the 
connectivity of the initially infected nodes. We used ini- 
tial conditions (i) with small connectivity k < 10 and 



(ii) with large connectivity k > 50 (the maximal connec- 
tivity of the network is 100 in this case). Figure ITTT b - ) 
shows that the condition (ii) is the main contribution to 
the small u < 1 values while it is the condition (i) which 
contributes mostly to the larger values u > 1. These re- 
sults can be easily explained as follows: We know from 
previous sections that the first steps of the epidemics pro- 
cess are concerning the hubs' infection. In the case (i), 
the path from the infected seeds to the hubs can be long 
and can thus induce a slower growth in i(t). In the op- 
posite case (ii), instead, the seeds are mainly located on 
the hubs since the beginning and the initial outbreak is 
generally faster that the average. This clear difference 
tells us that statistical deviations from the average time 
scale r can be traced back to the initial conditions of 
the epidemics. In this perspective, the specific network 
realizations appear to have a rather small influence on 
statistical fluctuations. This last feature is due to the in- 
trinsic large amount of connectivity fluctuations that each 
network realization possesses in the case of a very hetero- 
geneous topology. Additional fluctuations, such as those 
stemming from the specific realization and finite size, are 
thus a higher order correction to the intrinsic statistical 
fluctuations which are accounted for correctly in the the- 
oretical description. 
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Figure 11: Probability distribution of u — t' jr. (a) For a 
given network with random initial conditions (the BA network 
is obtained for N = 10 4 , m = 2, and 1,000 realizations), (b) 
Probability distribution of u for different initial condition: in 
full lines the nodes initially infected have a connectivity k < 10 
and in dashed lines k > 50 (here k max = 100). 



More generally, these results show that initial condi- 
tions induce fluctuations in the outbreak time which are 
never 'large' in the sense that P(u) decays rapidly, which 
implies that the time scale r is relevant even for a single 
realization of the network. This result is important since 
in real world situations, the disease spreads on one single 
network and it has to be checked whether results obtained 
on average are relevant or not. 

8. Conclusions 

In this paper we have provided a general picture of 
the effect of complex connectivity patterns in the spread- 
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ing dynamics of epidemic phenomena. We have shown 
how the large connectivity fluctuations present in a large 
class of population networks lead to a novel epidemio- 
logical framework in which both the epidemic threshold 
and the growth time-scale of the outbreak do not have 
an intrinsic value. Indeed, the virtual lack of epidemic 
threshold and instantaneous rising of the prevalence in in- 
finitely large networks corresponds to non-intrinsic quan- 
tities depending on the specific system size of the actual 
population. In addition, the dynamical structure of the 
spreading process is characterized by a hierarchical cas- 
cade from hubs to intermediate k and finally to small k 
vertices. That is, the infection first takes place on the 
subset of individuals with the largest number of contacts, 
and then progressively invades individuals with decreas- 
ing number of contacts. The emerging picture might be 
of practical importance in the implementation and as- 
sessment of dynamic control strategies. In particular, an 
efficient way to stop epidemics could rely on a dynami- 
cal deployment of containment measures that focuses on 
progressively changing classes of the population. More 
specifically, our results confirm the importance of control 
strategies targeting the hubs of the population but also 
highlight the fact that global surveillance is a major key 
aspect of epidemics control and that immunization strate- 
gies have to evolve with time during the different phases 
of the spread. 
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